#!/usr/bin/env python3
"""
糖尿病回归模型高级评估和医疗应用分析
===================================

本脚本提供糖尿病严重程度预测模型的深度评估，
包括医疗专用指标、模型可解释性、临床决策支持等分析。
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import metrics
from sklearn.inspection import permutation_importance
import joblib
import os
from scipy import stats
import warnings

warnings.filterwarnings('ignore')
plt.style.use('seaborn-v0_8')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

class DiabetesModelEvaluator:
    """糖尿病回归模型评估器"""
    
    def __init__(self, model_path=None, results_path=None):
        """
        初始化评估器
        
        Args:
            model_path: 模型文件路径
            results_path: 结果文件路径
        """
        self.model_path = model_path or '../results/diabetes_automl_model.pkl'
        self.results_path = results_path or '../results/evaluation_results.pkl'
        self.model = None
        self.results = None
        self.feature_names = [
            'age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6'
        ]
        self.feature_descriptions = {
            'age': '年龄',
            'sex': '性别', 
            'bmi': '体重指数',
            'bp': '平均血压',
            's1': '血清胆固醇',
            's2': '低密度脂蛋白',
            's3': '高密度脂蛋白',
            's4': '血清甘油三酯',
            's5': '血糖水平',
            's6': '血糖耐受指标'
        }
        
    def load_model_and_results(self):
        """加载模型和结果"""
        print("正在加载模型和评估结果...")
        
        if os.path.exists(self.model_path):
            self.model = joblib.load(self.model_path)
            print(f"✅ 模型加载成功: {self.model_path}")
        else:
            print(f"❌ 模型文件不存在: {self.model_path}")
            return False
            
        if os.path.exists(self.results_path):
            self.results = joblib.load(self.results_path)
            print(f"✅ 结果加载成功: {self.results_path}")
        else:
            print(f"❌ 结果文件不存在: {self.results_path}")
            return False
            
        return True
    
    def comprehensive_metrics_analysis(self):
        """全面的评估指标分析"""
        print("\\n📊 全面评估指标分析")
        print("=" * 60)
        
        test_metrics = self.results['test_metrics']
        train_metrics = self.results['train_metrics']
        
        # 创建指标对比表
        metrics_df = pd.DataFrame({
            '训练集': [
                train_metrics['r2_score'],
                train_metrics['rmse'],
                train_metrics['mae'],
                train_metrics['mape'],
                train_metrics['correlation'],
                train_metrics['clinical_acceptable_rate']
            ],
            '测试集': [
                test_metrics['r2_score'],
                test_metrics['rmse'],
                test_metrics['mae'],
                test_metrics['mape'],
                test_metrics['correlation'],
                test_metrics['clinical_acceptable_rate']
            ]
        }, index=['R²决定系数', 'RMSE', 'MAE', 'MAPE(%)', '相关系数', '临床可接受率(%)'])
        
        print("性能指标对比:")
        print(metrics_df.round(4))
        
        # 过拟合检测
        print(f"\\n🔍 过拟合检测:")
        r2_diff = train_metrics['r2_score'] - test_metrics['r2_score']
        if r2_diff > 0.1:
            print(f"⚠️  可能存在过拟合 (R²差值: {r2_diff:.3f})")
        elif r2_diff > 0.05:
            print(f"⚠️  轻微过拟合倾向 (R²差值: {r2_diff:.3f})")
        else:
            print(f"✅ 泛化性能良好 (R²差值: {r2_diff:.3f})")
        
        return metrics_df
    
    def residual_analysis(self):
        """残差分析"""
        print("\\n📈 残差分析")
        print("=" * 50)
        
        y_test = self.results.get('y_test', [])  # 需要在main.py中保存
        test_pred = self.results['test_predictions']
        
        if len(y_test) == 0:
            print("⚠️  无法进行残差分析 - 缺少真实测试标签")
            return
        
        residuals = y_test - test_pred
        
        # 残差统计
        print("残差统计信息:")
        print(f"- 残差均值: {np.mean(residuals):.4f}")
        print(f"- 残差标准差: {np.std(residuals):.4f}")
        print(f"- 残差中位数: {np.median(residuals):.4f}")
        print(f"- 残差范围: [{np.min(residuals):.2f}, {np.max(residuals):.2f}]")
        
        # 正态性检验
        shapiro_stat, shapiro_p = stats.shapiro(residuals)
        print(f"\\n残差正态性检验 (Shapiro-Wilk):")
        print(f"- 统计量: {shapiro_stat:.4f}")
        print(f"- p值: {shapiro_p:.4f}")
        if shapiro_p > 0.05:
            print("✅ 残差符合正态分布")
        else:
            print("⚠️  残差不符合正态分布")
        
        # 同方差性检验 (Breusch-Pagan test的简化版本)
        residuals_abs = np.abs(residuals)
        correlation_coef = np.corrcoef(test_pred, residuals_abs)[0, 1]
        print(f"\\n同方差性检验:")
        print(f"- 预测值与|残差|的相关系数: {correlation_coef:.4f}")
        if abs(correlation_coef) < 0.1:
            print("✅ 满足同方差性假设")
        else:
            print("⚠️  可能存在异方差性")
    
    def prediction_intervals_analysis(self):
        """预测区间分析"""
        print("\\n📏 预测区间分析")
        print("=" * 50)
        
        test_pred = self.results['test_predictions']
        test_metrics = self.results['test_metrics']
        
        # 基于RMSE的预测区间
        rmse = test_metrics['rmse']
        
        print("预测区间设置:")
        print(f"- 68%置信区间: ±{rmse:.2f}")
        print(f"- 95%置信区间: ±{1.96 * rmse:.2f}")
        print(f"- 99%置信区间: ±{2.58 * rmse:.2f}")
        
        # 实际覆盖率分析（如果有真实值）
        y_test = self.results.get('y_test', [])
        if len(y_test) > 0:
            residuals = np.abs(y_test - test_pred)
            
            coverage_68 = np.mean(residuals <= rmse) * 100
            coverage_95 = np.mean(residuals <= 1.96 * rmse) * 100
            coverage_99 = np.mean(residuals <= 2.58 * rmse) * 100
            
            print(f"\\n实际覆盖率:")
            print(f"- 68%区间覆盖率: {coverage_68:.1f}%")
            print(f"- 95%区间覆盖率: {coverage_95:.1f}%")
            print(f"- 99%区间覆盖率: {coverage_99:.1f}%")
    
    def clinical_performance_analysis(self):
        """临床性能分析"""
        print("\\n🏥 临床性能分析")
        print("=" * 50)
        
        test_metrics = self.results['test_metrics']
        y_test = self.results.get('y_test', [])
        test_pred = self.results['test_predictions']
        
        # 临床风险分层分析
        if len(y_test) > 0:
            # 定义风险等级
            risk_thresholds = [np.percentile(y_test, 33), np.percentile(y_test, 67)]
            
            def get_risk_level(value):
                if value <= risk_thresholds[0]:
                    return '低风险'
                elif value <= risk_thresholds[1]:
                    return '中风险'
                else:
                    return '高风险'
            
            true_risk = [get_risk_level(v) for v in y_test]
            pred_risk = [get_risk_level(v) for v in test_pred]
            
            # 风险分层准确性
            risk_accuracy = np.mean([t == p for t, p in zip(true_risk, pred_risk)]) * 100
            print(f"风险分层准确性: {risk_accuracy:.1f}%")
            
            # 风险等级混淆矩阵
            from sklearn.metrics import confusion_matrix
            risk_labels = ['低风险', '中风险', '高风险']
            cm = confusion_matrix(true_risk, pred_risk, labels=risk_labels)
            
            print(f"\\n风险分层混淆矩阵:")
            risk_cm_df = pd.DataFrame(cm, index=risk_labels, columns=risk_labels)
            print(risk_cm_df)
        
        # 临床决策阈值分析
        print(f"\\n📋 临床决策建议:")
        
        clinical_rate = test_metrics['clinical_acceptable_rate']
        if clinical_rate >= 85:
            decision = "推荐临床应用"
            confidence = "高置信度"
        elif clinical_rate >= 70:
            decision = "谨慎临床应用"
            confidence = "中等置信度"
        else:
            decision = "不推荐临床应用"
            confidence = "低置信度"
        
        print(f"- 决策建议: {decision}")
        print(f"- 置信度级别: {confidence}")
        print(f"- 临床可接受率: {clinical_rate:.1f}%")
        
        # 高风险患者预测性能
        if len(y_test) > 0:
            high_risk_threshold = np.percentile(y_test, 75)
            high_risk_true = y_test >= high_risk_threshold
            high_risk_pred = test_pred >= high_risk_threshold
            
            if np.sum(high_risk_true) > 0:
                sensitivity = np.sum(high_risk_true & high_risk_pred) / np.sum(high_risk_true)
                specificity = np.sum(~high_risk_true & ~high_risk_pred) / np.sum(~high_risk_true)
                
                print(f"\\n高风险患者识别能力:")
                print(f"- 敏感性 (召回率): {sensitivity:.3f}")
                print(f"- 特异性: {specificity:.3f}")
    
    def model_interpretability_analysis(self):
        """模型可解释性分析"""
        print("\\n🔬 模型可解释性分析")
        print("=" * 50)
        
        # Auto-sklearn模型构成分析
        print("模型集成构成:")
        try:
            models_info = str(self.model.show_models())
            print("模型详情:")
            print(models_info[:500] + "..." if len(models_info) > 500 else models_info)
        except Exception as e:
            print(f"无法获取模型详情: {e}")
        
        # 特征重要性分析（如果可能）
        print(f"\\n特征重要性估计:")
        print("注: Auto-sklearn是黑箱集成模型，特征重要性需要通过其他方法估计")
        
        # 基于相关性的特征重要性估计
        test_metrics = self.results['test_metrics']
        if 'feature_correlations' in self.results:
            correlations = self.results['feature_correlations']
            print("基于相关性的特征重要性排序:")
            sorted_features = sorted(correlations.items(), key=lambda x: abs(x[1]), reverse=True)
            for i, (feature, corr) in enumerate(sorted_features[:5], 1):
                desc = self.feature_descriptions.get(feature, feature)
                print(f"{i}. {desc}: {abs(corr):.3f}")
    
    def comparative_analysis(self):
        """对比分析"""
        print("\\n⚖️  对比分析")
        print("=" * 50)
        
        test_metrics = self.results['test_metrics']
        
        # 与基准模型比较
        print("与基准模型性能比较:")
        print("(基准: 线性回归、随机预测)")
        
        # 简单基准 - 均值预测
        y_test = self.results.get('y_test', [])
        if len(y_test) > 0:
            mean_baseline_r2 = 0.0  # 均值预测的R²为0
            current_r2 = test_metrics['r2_score']
            improvement = current_r2 - mean_baseline_r2
            
            print(f"- 相对均值预测提升: {improvement:.3f}")
            print(f"- 性能提升倍数: {current_r2 / 0.001 if current_r2 > 0 else 'N/A':.1f}x")
        
        # 医疗标准比较
        print(f"\\n医疗AI标准评估:")
        r2_score = test_metrics['r2_score']
        
        if r2_score >= 0.8:
            grade = "A级 (优秀)"
        elif r2_score >= 0.6:
            grade = "B级 (良好)"
        elif r2_score >= 0.4:
            grade = "C级 (可接受)"
        else:
            grade = "D级 (需改进)"
        
        print(f"- 模型等级: {grade}")
        print(f"- R²分数: {r2_score:.3f}")
    
    def create_advanced_visualizations(self):
        """创建高级可视化图表"""
        print("\\n正在生成高级评估图表...")
        
        y_test = self.results.get('y_test', [])
        test_pred = self.results['test_predictions']
        
        if len(y_test) == 0:
            print("⚠️  缺少测试标签，跳过部分可视化")
            return
        
        # 创建综合评估图表
        fig, axes = plt.subplots(2, 3, figsize=(18, 12))
        fig.suptitle('糖尿病回归模型深度评估分析', fontsize=16, fontweight='bold')
        
        # 1. 预测准确性散点图
        ax1 = axes[0, 0]
        ax1.scatter(y_test, test_pred, alpha=0.6, color='steelblue')
        ax1.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
        ax1.set_xlabel('真实值')
        ax1.set_ylabel('预测值')
        ax1.set_title('预测准确性')
        
        # 添加统计信息
        r2 = self.results['test_metrics']['r2_score']
        rmse = self.results['test_metrics']['rmse']
        ax1.text(0.05, 0.95, f'R² = {r2:.3f}\\nRMSE = {rmse:.2f}', 
                transform=ax1.transAxes, bbox=dict(boxstyle="round", facecolor='wheat'))
        
        # 2. 残差分布
        ax2 = axes[0, 1]
        residuals = y_test - test_pred
        ax2.hist(residuals, bins=20, alpha=0.7, color='orange', edgecolor='black')
        ax2.set_xlabel('残差')
        ax2.set_ylabel('频次')
        ax2.set_title('残差分布')
        ax2.axvline(x=0, color='r', linestyle='--')
        
        # 3. 预测区间图
        ax3 = axes[0, 2]
        sorted_indices = np.argsort(test_pred)
        sorted_pred = test_pred[sorted_indices]
        sorted_true = y_test[sorted_indices]
        rmse = self.results['test_metrics']['rmse']
        
        ax3.plot(sorted_pred, sorted_true, 'o', alpha=0.6, label='实际值')
        ax3.plot(sorted_pred, sorted_pred, 'r-', label='完美预测')
        ax3.fill_between(sorted_pred, sorted_pred - 1.96*rmse, sorted_pred + 1.96*rmse, 
                        alpha=0.2, color='gray', label='95%预测区间')
        ax3.set_xlabel('预测值')
        ax3.set_ylabel('真实值')
        ax3.set_title('预测区间分析')
        ax3.legend()
        
        # 4. 绝对误差分布
        ax4 = axes[1, 0]
        abs_errors = np.abs(residuals)
        ax4.boxplot(abs_errors)
        ax4.set_ylabel('绝对误差')
        ax4.set_title('绝对误差分布')
        ax4.grid(True, alpha=0.3)
        
        # 5. 相对误差分析
        ax5 = axes[1, 1]
        relative_errors = np.abs((y_test - test_pred) / y_test) * 100
        ax5.hist(relative_errors, bins=20, alpha=0.7, color='purple', edgecolor='black')
        ax5.set_xlabel('相对误差 (%)')
        ax5.set_ylabel('频次')
        ax5.set_title('相对误差分布')
        ax5.axvline(x=20, color='r', linestyle='--', label='临床阈值')
        ax5.legend()
        
        # 6. 性能趋势分析
        ax6 = axes[1, 2]
        # 按预测值分组分析性能
        pred_bins = np.linspace(test_pred.min(), test_pred.max(), 10)
        bin_centers = []
        bin_mae = []
        
        for i in range(len(pred_bins)-1):
            mask = (test_pred >= pred_bins[i]) & (test_pred < pred_bins[i+1])
            if np.sum(mask) > 0:
                bin_centers.append((pred_bins[i] + pred_bins[i+1]) / 2)
                bin_mae.append(np.mean(np.abs(y_test[mask] - test_pred[mask])))
        
        ax6.plot(bin_centers, bin_mae, 'o-', color='green')
        ax6.set_xlabel('预测值区间中心')
        ax6.set_ylabel('平均绝对误差')
        ax6.set_title('不同预测值范围的性能')
        ax6.grid(True, alpha=0.3)
        
        plt.tight_layout()
        
        # 保存图表
        os.makedirs('../results', exist_ok=True)
        plt.savefig('../results/advanced_evaluation.png', dpi=300, bbox_inches='tight')
        
        print("✅ 高级评估图表已保存至 ../results/advanced_evaluation.png")
    
    def generate_clinical_deployment_report(self):
        """生成临床部署报告"""
        print("\\n🏥 临床部署建议报告")
        print("=" * 60)
        
        test_metrics = self.results['test_metrics']
        
        print("📋 部署准备状态评估:")
        
        # 性能检查清单
        checklist = {
            '模型准确性': test_metrics['r2_score'] >= 0.5,
            '预测稳定性': test_metrics['rmse'] <= 40,
            '临床可接受性': test_metrics['clinical_acceptable_rate'] >= 70,
            '相关性强度': test_metrics['correlation'] >= 0.6
        }
        
        passed_checks = sum(checklist.values())
        total_checks = len(checklist)
        
        for check, status in checklist.items():
            status_symbol = "✅" if status else "❌"
            print(f"{status_symbol} {check}")
        
        print(f"\\n📊 整体就绪度: {passed_checks}/{total_checks} ({passed_checks/total_checks*100:.0f}%)")
        
        # 部署建议
        if passed_checks >= 3:
            print("\\n🚀 部署建议: 推荐部署")
            print("适用场景:")
            print("- 辅助诊断工具")
            print("- 患者风险分层")
            print("- 治疗效果预估")
        else:
            print("\\n⚠️  部署建议: 需要优化")
            print("改进方向:")
            print("- 收集更多高质量数据")
            print("- 优化特征工程")
            print("- 调整模型超参数")
        
        # 监控指标建议
        print("\\n📊 部署后监控指标:")
        print("- 预测准确性监控 (每周)")
        print("- 数据漂移检测 (每月)")
        print("- 临床反馈收集 (持续)")
        print("- 模型性能评估 (每季度)")
        
        # 风险评估
        print("\\n⚠️  风险评估和缓解:")
        print("潜在风险:")
        print("- 数据分布变化影响预测准确性")
        print("- 极端值预测可能不准确")
        print("- 模型黑箱性限制可解释性")
        
        print("\\n缓解措施:")
        print("- 建立数据质量监控系统")
        print("- 设置预测值合理性检查")
        print("- 提供预测置信度信息")
        print("- 定期模型重训练和验证")

def main():
    """主函数"""
    print("🔬 糖尿病回归模型深度评估")
    print("=" * 50)
    
    # 创建评估器
    evaluator = DiabetesModelEvaluator()
    
    try:
        # 加载模型和结果
        if not evaluator.load_model_and_results():
            print("❌ 无法加载必要文件，请先运行 main.py")
            return
        
        # 执行全面评估
        evaluator.comprehensive_metrics_analysis()
        evaluator.residual_analysis()
        evaluator.prediction_intervals_analysis()
        evaluator.clinical_performance_analysis()
        evaluator.model_interpretability_analysis()
        evaluator.comparative_analysis()
        evaluator.create_advanced_visualizations()
        evaluator.generate_clinical_deployment_report()
        
        print("\\n🎉 模型深度评估完成！")
        print("📂 评估图表已保存至 ../results/ 目录")
        
    except Exception as e:
        print(f"❌ 评估过程出错: {str(e)}")
        raise

if __name__ == "__main__":
    main()